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SUMMARY 


The  paper  addresses  the  finite  element  method  with  a-posteriori  error 
estimation  for  elements  of  degree  p  -  1  and  p  -  2.  It  gives  the  formu¬ 
lae  for  the  error  indicators  and  error  estimators.  Basic  mathematical 
characterization  of  the  estimators  are  given  and  it  is  shown  that  the 
estimators  for  p  -  1  and  p  -  2  have  different  structures.  Numerical 
examples  show  the  effectivity  of  the  approach  and  the  high  quality  of  the 


estimator . 


1 .  INTRODUCTION 


During  the  last  few  years  a  significant  progress  has  been  achieved 
in  the  theory  and  the  implementation  of  the  adaptive  procedures  and 
a-posteriori  error  estimation  for  the  finite  element  method.  For  the 
survey  of  the  today's  state  of  the  art  of  the  quality  assessment  of  finite 
element  solutions  we  refer  to  [1].  Adaptive  finite  element  codes  FEARS 
[2],  PLTGM  [3].  EXPDES  [4]  were  developed.  For  the  theoretical  and 
implementational  aspects  related  to  the  mentioned  codes  we  refer  to  [5]  to 
[10].  Various  approaches  for  obtaining  a-posteriori  error  estimates  were 
recently  developed.  See,  for  example,  [1]  -  [15]  and  other  references 
mentioned  in  [1]. 

The  codes  and  the  theory  mentioned  above  are  for  the  finite  element 
method  with  elements  of  degree  one.  In  this  paper  we  are  addressing  the 
questions  arising  from  the  use  of  elements  of  degree  two.  We  will  concen¬ 
trate  on  the  comparison  between  the  approaches  based  on  the  elements  of 


degree  one  and  two. 


2.  THE  MODEL  PROBLEM 


Let  us  consider  the  model  problem  of  the  Laplace  equation 

(2.1.)  L(u)  -  -Au  -  f  on  Q  c 

with  the  boundary  conditions 

0 

(2.2a)  u  -  0  on  3  8 

(2.2b)  -  t  on  31 n 

on 

where  n  is  the  outward  normal  to  3D. 

The  solution  uQ  of  the  problem  (2.1)  (2.2)  is  understood  in  the 
standard  weak  sense:  Uq  €  Hq(8)  -  {u|u  €  H^(Q),  u  -  0  on  3°8}  such 
that 

(2.3)  B(uq,v)  -  j  (VuQ*Vv)dx  -  /  fvdx  +  J  tvds 

Q  B  31  8 

holds  for  all  v  €  Hq(8).  By  ^(8)  we  denote  the  standard  Sobolev 
2 

space.  By  |u|E  -  B(u,u)  we  denote  the  energy  norm  of  function 
u  €  Hj(Q). 

Remark.  Although  we  restricted  ourself  to  a  very  special  case  (2.1), 
(2.2)  the  theory  and  conclusions  hold  for  the  general  elliptic  equation  on 
bounded  domains  in  R  . 

We  will  consider  two  illustrative  model  problems. 

Problem  1 .  The  domain  8  and  the  boundary  conditions  are  shown  in 


Fig.  2.1. 


Fig.  2.1.  Scheme  of  Problem  1. 


3°Q  consists  of  OA  and  on  3^  t  is  defined  so  that  the  exact 
solution  uq  is 

uQ  -  0.070075 sin  0/2 

and 

|u0|E  -  0.657878  *  10*1. 

By  (r,0)  we  denote  the  polar  coordinates. 

Problem  2.  The  (slit)  domain  0  and  the  boundary  conditions  are 
shown  in  Fig.  2.2. 

3°n  consists  of  upper  side  0A  (marked)  and  on  t  is  defined 

so  that  the  exact  solution  Uq  is 

u0  -  0.0635659rl/*  sin  0/4 

|u0|E  -  0.579^28  *  10-1 

Problem  1  and  2  have  obviously  solution  with  different  strength  of 


singularity. 
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3.  THE  FINITE  ELEMENT  METHOD 


We  will  consider  the  finite  element  method  with  the  bilinear  and 
biquadratic  elements  on  the  meshes  which  are  used  in  the  Program  FEARS. 
These  meshes  have  regular  and  irregular  nodel  points.  The  values  of  the 
finite  element  solution  in  the  irregular  nodal  points  is  determined  by  the 
requirement  that  the  elements  are  conforming.  Fig.  3-1a,b  show  the  meshes 
with  20  elements  for  the  model  problem  1  and  bilinear  and  biquadratic  ele¬ 
ments.  The  regular  nodel  points  are  marked  by  •,  while  the  irregular 
points  are  marked  by  *.  We  will  not  introduce  here  a  formal  definition 
and  refer  to  [2],  [6],  [7]. 

We  have  shown  the  meshes  only  for  the  special  case.  They  are  in 
general  defined  for  domains  with  curved  boundaries.  See  [2],  [6]. 

The  finite  element  solution  upE  is  defined  in  the  standard  way. 

By  e  we  denote  the  error  of  the  finite  element  solution  uFE,  e  - 

uFE  “  u* 


< 

4 

t 

\ 

■ - * 

c 

Fig.  3.1a.  The  regular  and  irregular  nodal  points  for  p  -  1 . 


nwati 
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4.  THE  ADAPTIVE  FINITE  ELEMENT  METHOD  A-POSTERIORI  ERROR 
ESTIMATION  AND  THE  QUALITY  ASSESSMENT 

We  will  distinguish  between  the  feedback  and  adaptive  finite  element 
method.  A  feedback  method  utilizes  the  currently  available  information 
for  steering  the  computational  process.  The  adaptive  approach  is  a  feed¬ 
back  approach  which  has  clearly  defined  optimal  properties.  Hence  the 
assessment  of  the  quality  of  a  feedback  approach  is  relative  to  the 
criterion  of  the  optimality.  We  will  not  go  here  into  details  but  refer 
to  [16],  [17]. 

We  will  say  that  the  feedback  is  adaptive  with  respect  to  the 
convergence  rate  measure  if  the  rate  is  the  same  as  the  theoretically 
best  one  among  all  the  meshes  with  the  same  number  of  degrees  of  freedom 
N  and  the  same  degrees  of  elements  p. 

For  our  model  problems  the  feedback  approach  will  be  adaptive  if  the 
(asymptotic)  rage  is  0(N  2  )  for  the  elements  of  degree  p  »  1  and  is 
0(N-1)  for  the  elements  of  degree  p  =  2  in  both  cases,  i.e.  the  rate  is 
independent  of  the  strength  of  the  singularity.  We  note  that  the  same 
rate  is  obtained  when  the  solution  is  smooth. 

The  quality  of  the  a-posteriori  error  estimator  e  is  measured  by 
its  efx'ecti vity  index  Q  =  e/|e|g  where  [e|E  is  energy  norm  of  the 
error  of  the  finite  element  solution.  In  practice  we  have  to  require 
that  1 0-1  |  <,  .2  (say)  when  the  accuracy  of  Upg  is  in  the  range  of  the 
engineering  requirement  (i  10$  say).  In  addition  we  prefer  that  0  >  1 
so  that  the  true  error  is  overestimated  rather  than  underestimated.  This 
is  important  because  the  error  estimator  is  used  as  the  stopping 
criterion.  Obviously  if  6  >  2  or  9  <  \  ,  then  the  error  estimator  is 
practically  unacceptable. 
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5.  THE  ERROR  INDICATORS  AND  ESTIMATORS 

The  error  estimator  e  is  computed  through  the  error  indicators 
n(A)  of  the  single  element  A.  Then  the  error  estimator  (which 
approximate  |e|E)  is  given  by 

(5.1)  e2  -  l  n( A)2 

A 

where  the  sum  is  taken  over  all  elements. 

The  error  indicators  are  also  used  for  the  steering  the  feedback 
finite  element  method.  We  refer  to  [18]  for  the  detailed  analyses  of 
various  steering  strategies.  The  basic  strategy  is  to  refine  all  elements 
having  the  error  indicators  over  the  threshold  t,  for  example,  t  * 

Y  max  n( A)  where  0  <  Y  ^  1  is  a-priori  chosen  and  the  maximum  is  taken 
over  all  elements. 

The  code  FEARS  uses  basically  this  type  of  strategy  but  more 
sophisticated  to  minimize  computer  time  by  avoiding  repeated 
recomputation.  The  experience  with  FEARS  shows  that  the  total  computer 
time  is  2-3  times  the  computer  time  solving  the  final  mesh  only. 

We  note  that  error  indicators  which  lead  to  completely  unacceptable 
quality  of  the  error  estimators  can  still  perform  relatively  well  in 
steering  the  feedback  approach. 

Let  us  now  define  the  error  indicators  for  the  bilinear  and 
biquadratic  elements.  For  simplicity  of  the  exposition,  we  will  consider 
only  an  uniform  patch  of  the  mesh.  The  data  management  of  the  adaptive 
solver  is  based  on  the  tree  structure  (see  [7])  and  mesh  is  defined  in  a 
recursive  way.  Hence,  every  element  has  a  "father"  which  has  four 


"sons."  Fig.  5.1  shows  uniform  mesh  with  element  A.j',  i  -  1,2, 3, 4 
belonging  to  one  "father"  (A0)  shown  also  in  the  figure. 


Fig.  5.1.  The  element  and  his  "father". 

A.  .  Error  Indicator  for  the  bilinear  elements. 

We  will  define  the  error  indicator  as  used  in  FEARS  but  simplified 
for  the  special  problem  under  consideration. 

Let 

2 

(5.2a)  n!(A?)  »  |A|2  l  (  max  r  (A.))2 

J-1  i-1 . 4  j  1 

(5.2b)  -  | A|s  |b|2  0 

L2(  V 

where  Ja|  is  the  length  of  the  element  side,  rx^(Aj)  is  the  jump  of 
the  derivative  in  the  x,  direction  in  the  vertex  A4  of  the  element 


(5.3a)  R  -  f  +  AupE  -  p 

where  p  is  such  that  J  R  dx  -  0. 

Ao 

Now  we  define  i 

(5.4a)  n2(A?)  -  n?(A?)  ♦  n|( A?) . 

It  can  be  proven  (see  [6])  that  ^(A®)  <<  hi  (A®)  and  hence  it  can  be 
neglected  and  define 

(5.4b)  n2(  A°)  -  nf(Af). 

Experience  with  FEARS  shows  that  in  fact  the  use  of  n  and  n 
practically  does  not  change  the  estimator  in  [6]. 

Remark  1 .  If  the  mesh  is  uniform  and  the  solution  is  quadratic,  then 
the  error  estimator  is  exact. 

Remark  2.  We  could  use  other  equivalent  estimators,  for  example, 

O 

replace  by  the  integral  of  the  square  of  the  jump  of  the  derivative 

of  the  finite  element  solution  over  the  boundary  of  A®  (multiplied  by 
|a|).  Obviously  n2(A^)  is  the  value  of  this  integral  computed  by  a 
quadrative.  Nevertheless,  experience  shows  that  for  small  accuracy  (large 
error)  the  estimator  based  on  the  integral  tends  to  underestimate  the 
error  and  hence  the  definition  of  the  indicator  n1  corrects  this 
undesirable  and  enlarges  the  range  of  the  asymptotic  validity  of  the 
situation  of  the  estimator  (see  al3o  numerical  results  in  section  7). 


.A 


11 


aJ  WJ  , 


B.  Error  Indicator  for  biquadratic  elements 


Let  us  assume  that  the  center  of  the  element  lies  in  the  origin. 


Define 


(5.5a) 


GUj) 


(5.5b) 


(5.6) 


>i  -  *?  -  if-- 


V1  -  xlM2 


/  (f  vt  -  TUpg  •W1)dx, 


GtVi), 


n?u?>  ■  (0(v‘,>: 

Ml  1-1 


i  -  1,2 


i  “1,2 


(The  integrals  in  (5.5a)  are  numerically  computed.) 


(5.7) 


nf(A°)  -  |A|  l  J(rx  )2ds 

j-1  j 


where  the  integral  is  taken  over  the  sides  of  the  elements  perpendicular 

to  x,  and  rv  are  as  before  the  jump  in  the  derivatives. 

J  Xj 


(5.8) 


where 


-  |4|2  |R|2  o 

L2(V 


R  “  f  +  AU 


U  *  UFE  +  b1x161  +  b2x262 


and  are  given  by  (5.5b).  Then  we  define 

n2(  A?)  -  I  nf(  A?) 

i-1 

and  e  is  defined  by  (5.1). 

It  can  be  shown,  see  [19],  that  if  the  solution  is  sufficiently 
smooth,  then  n2  311(1  ^3  can  136  neglected  with  respect  to  ni  and  hence 
we  can  use  the  indicator 

n(A°)  -  n^A®) 

as  in  the  case  for  p  -  1.  Nevertheless,  we  see  a  significant  difference, 
namely  that  in  the  case  p  -  1  the  principal  part  of  the  indicator  is  the 

"boundary  line  part"  (jumps),  while  for  p  -  2  it  is  the  "volume"  part. 

From  the  implementational  point  of  view,  the  "volume"  indicator  is  more 
advantageous . 

In  the  case  of  the  smooth  solution  uQ  and  reasonable  meshes  the 
indicators  for  p  -  1  and  p  -  2  are  of  the  same  magnitude  and  the 
estimators  are  of  the  high  quality,  namely 

|e|£  -  e(1  +  0(e) ) , 

where  the  term  0(e)  is  very  small.  Nevertheless,  if  the  solution  is 
unsmooth  as  in  our  model  problems,  then  the  indicators  are  of  very  dif¬ 
ferent  magnitude  when  the  elements  are  adjacent  to  the  origin  where  the 
solution  is  singular  and  the  error  is  in  these  elements  strongly  under¬ 
estimated.  This  is  especially  important  in  the  case  p  -  2  because  the 

error  indicators  in  the  area  where  the  solution  is  smooth  are  much  smaller 

than  for  p  -  1.  This  makes  the  contribution  of  the  largest  indicators  to 
the  error  estimator  much  larger. 
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The  elements  which  are  in  the  area  of  the  singularity  of  the  solution 
are  refined  during  the  adaptive  mesh  construction  and  are  located  in  the 
terminal  points  of  the  tree  data  structure.  These  elements  are  quadruples 
of  the  elements  belonging  to  one  "father".  Fig.  5.2  shows  in  the  left 
corner  such  a  quadruple  and  its  "father"  and  "grandfather". 


Fig.  5.2.  Scheme  of  the  terminal  element. 

Hence,  we  will  modify  the  error  indicators  in  these  terminal  quadruples  in 
dependence  on  the  distribution  of  the  error  indicators.  (This  modifica¬ 
tion  can  be  implemented  in  a  very  effective  way.) 

A 

The  modified  estimators  n  are  constructed  from  n  as  follows:  the 
four  indicators  of  the  quadruple  of  the  elements  are  ordered  by  the 
magnitude.  For  the  sake  of  the  explanation,  assume  that  n2(A®)  £  ^(A^) 

>  n2( Ag)  >  n2(Ajj)  (see  Fig.  5.1). 

If  n2(A®)  >  Yn2(Aj)  with  Y  -  8,  then 

n2(  A^)  -  n2(A?)F(R) 

n2(A°)  -  n2(A°),  i  -  2,3,H 


F(R) 


0.031  605R  ♦  1 


for  0  i  R  £  23*1 


-  3. 

48l3e3*762* 

10  3r  for  234 

£  R  £  1000 

-  150 

for  R  > 

1000 

with  R 

2 ,  .0. 
n  A1 

-  -5 — fr-  where  An  is 
m2(a2) 

the  diagonal  element  to  A® 

Yn2(  A®) 

and  n2(A°) 

>  Bn2(A§) 

with  $  -  100, 

then 

n2(A?)  - 

GfRjJr^U?) 

n2(A0)  - 

G(R3)n2(A^) 

n2( A°)  - 

n2( a|) 

n2ug)  - 

n2( A®) 

If  n2(A?) 


where 


and 


0.003499R  ♦  1 

for 

R  i  700 

1.7132eR’10  3 

for 

700  i  R  £  4000 

93-54 

for 

4000  i  R 

2,  .0. 
n  ( A1 ) 

D  m 

2.  .0. 
n  ( A^) 

2,  .0/ 
n  (a4) 

R3 

2.  A0. 
n  (a2) 

(where  as  before  In  the  denominator  is  the  indicator  of  the  diagonally 
opposite  element).  If  n2(A3)  <  Sh^CA^).  then 


The  above  formulae  are  characterizing  the  singularity  of  a  point  character 


( F( R) )  (e.g.  of  the  type  occurring  in  our  model  problems)  or  the  line 
singularity  along  the  element  side  (G(R)).  The  formulae  were  developed  by 
analyzing  theoretically  and  numerically  these  cases  and  their  validity  was 
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6.  SOME  BASIC  PROPERTIES  OF  THE  ESTIMATORS 

Let  us  mention  the  essential  properties  of  the  estimator.  We  list 
the  basic  assumptions  but  not  completely  accurately  because  their  com¬ 
plexity.  We  refer  to  [6]  and  [19]  to  precise  formulations. 

a)  The  mesh  is  a  K-mesh.  Roughly  speaking,  the  mesh  is  a  K-mesh  if 
there  are  not  too  many  irregular  nodal  points  between  the  regular  ones. 

b)  The  mesh  is  nearly  equilibrated.  This  means  that  the  error 
indicators  are  nearly  equal  (except  possibly  small  number  of  exceptions). 

c)  |e|g  2  ChP  where  h  is  the  maximal  element  size.  The 
conditions  (a,b)  usually  lead  to  c. 

d)  The  meshes  are  patchwise  uniform.  This  roughly  means  that  large 
pieces  of  meshes  have  uniform  character. 

e)  The  mesh  is  such  that  (roughly  speaking)  the  major  part  of  the 
error  is  in  the  areas  where  the  solution  has  no  singularity  (due  to  proper 
refinement)  (usually  b  leads  to  this  assumption). 

Under  the  above  assumptions,  the  following  theorem  holds  (see  [19]). 

Theorem  6.1.  For  any  t  >  0  and  h  £  )1q(t),  t  sufficiently  small 

1 0  -  1  |  <;  Mx 

where  the  constant  M  is  independent  of  x. 

Theorem  6.1  is  valid  for  p  *  1  and  p  -  2.  Also  for  not  modified 
estimator.  Nevertheless,  Mt  could  be  small  only  for  very  small  h, 
i.e.  for  very  high  accuracy  which  is  not  practically  relevant.  Hence  the 
modified  estimator  is,  as  we  will  see,  essential  for  practical  computa¬ 


tion  . 


7-  NUMERICAL  EXAMPLES 

We  will  illustrate  the  performance  of  the  feedback  finite  element 
method  and  the  error  estimator  based  on  the  error  indicators  introduced  in 
previous  sections.  Let  us  summarize  typical  behavior  one  can  observe 

1)  Fig.  7.1  shows  the  values  of  the  modified  error  Indicators  n  (a) 
for  Problem  1 ,  biquadratic  elements  (p=2)  and  various  meshes  with  M 
elements  and  N  degrees  of  freedom. 

o 

All  indicators  are  of  the  same  magnitude  (10  )  except  the  ones 

adjoined  to  the  singularity  point  which  are  of  much  larger  magnitude 
(10“5). 

The  indicators  of  the  elements  which  are  not  refined  do  not  change 
significantly  by  the  refinement  of  other  elements. 

2)  The  maximal  error  indicator  is  decreased  by  the  factor  j  for 
problem  1  (by  /2  /2  for  problem  2)  when  the  element  is  refined  by 

dividing  it  into  four  elements.  This  factor  is  directly  related  to  the 

y  y 

strength  of  the  singularity  (r2,  r'1*)  of  the  solution.  (Note  that  this 
factor  i3  independent  of  p).  For  p  «  2  and  problem  1  the  total  error 
e  is  essentially  the  maximal  error  indicator  until  the  indicator  is  of 

_7 

the  order  10  which  occurs  when  the  smallest  element  is  of  order 
2_1®,  i.e.  of  the  order  - and  the  error  |e|E  is  of  the  magnitude 
1J  of  |UIE*  Hence  the  rate  is  exponential.  Fig.  7.2a,  respectively 
7.2b,  shows  in  the  semilogarithmic  resp.  double  logarithmical  scale  the 
behavior  of  the  error  |e|E-  Similar  behavior  can  be  observed  in  the  case 
of  problem  2.  See  Fig.  7.3a,b.  In  Fig.  7.2b  and  7.3b  the  slope  of  the 
maximal  asymptotic  rate  N-''2  and  N-1  for  the  elements  of  degree 
p  -  1  and  p  -  2  is  shown.  This  rate  is  achieved  for  adaptively 
constructed  meshes  for  high  accuracies  independently  of  the  strength  of 
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the  singularity.  We  see  that  for  p  =  1  this  rate  is  achieved  in  the 
range  of  engineering  accuracy  while  for  p  -  2  the  rate  is  in  this  range 
still  exponential  and  the  asymptotic  rate  N-1  is  achieved  for  much 
higher  accuracies. 

3)  When  the  indicator  of  the  element  adjoined  to  the  origin  will  be 
of  order  of  the  other  indicators,  then  many  elements  will  be  refined.  The 

p 

indicators  n  (A)  in  the  area  where  the  solution  is  smooth  will  decrease 
by  the  magnitude  2^  -  64  for  p-  2  and  2^  »  16  for  p  -  1.  This 
refinement  is  then  followed  by  a  series  of  refinement  of  the  elements 
adjoined  to  the  singularity.  Fig.  7.2  shows  clearly  this  character  for  p 
-  1.  For  p  -  2  the  refinement  of  many  elements  did  not  occur  in  the 
figure . 

4)  The  refinement  by  the  bisection  procedure  is  not  an  optimal  one 
for  the  convergence.  Nevertheless,  it  is  very  advantageous  for  the  data 
treatment.  We  refer  to  [20]  for  the  analyses  of  the  optimal  refinement  in 
one  dimension. 

5)  An  effective  feedback  procedure  has  (especially  for  p  »  2) 
anticipated  the  character  of  the  proper  refinement  so  that  the  proper  mesh 
is  designed  quickly.  It  is  obvious  that  the  refinement  character  can  be 
recognized.  In  the  program  FEARS  this  procedure  is  called  a  short  pass 
procedure . 

6)  The  quality  of  the  error  estimator  is  high.  Fig.  7.2  shows  the 
error  and  the  error  indicator  for  p  =  2.  In  Table  7.1  we  show  typical  data 
for  the  meshes  with  8-26  elements  for  Problem  1.  We  clearly  see  that  the 
effectivity  index  for  the  nonmodified  estimator  is  unacceptable  although  it 
increases  with  |e|E  -*■  0  (and  theoretically  converges  to  1 ) .  On  the  other 


hand,  the  modified  estimator  has  effectivity  index  of  high  quality. 


7)  Comparing  performance  of  the  bilinear  (p  ■  1 )  and  biquadratic 
mesh,  we  see  that  in  the  presence  of  a  singularity  the  bilinear  elements 


perform  better  for  low  accuracy  while  higher  accuracy  cannot  be 
practically  achieved  by  bilinear  elements.  If  the  solution  is  smooth, 
then  the  biquadratic  elements  perform  better  in  larger  range  of  the 
accuracy  (see  also  Fig.  7.3a,b).  The  performance  of  biquadratic  elements 
could  be  improved  if  stronger  refinement  than  based  on  bisection  could  be 
effectively  made. 
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Fig.  7.1a.  Values  of  the  modified  error  indicators  n(A)  for  M  -  8,14. 
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Fig.  7.2a.  Accuracy  of  the  finite  element  solution  of  Problem  1  for  the 
element  of  degree  p  »  1,2  (semilog  scale). 
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Fig.  7.2b.  Accuracy  of  the  finite  element  solution  of  Problem  1  for  the 
elements  of  degree  p  -  1,2  (log  log  scale). 


Fig.  7.3a.  Accuracy  of  the  finite  element  solution  of  Problem  2  for  the 
elements  of  degree  p  -1,2  (semilog  scale). 
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8.  CONCLUSIONS 

We  summarize  now  briefly  our  conclusions. 

a)  The  error  indicators  have  different  structures  for  even  and  odd 
degree  of  elements.  The  even  degrees  are  preferable  from  the  implemen- 
tational  aspects. 

b)  The  range  of  asymptotic  exactness  of  the  estimator  has  to  be  in¬ 
creased  by  the  modification  of  the  indicators  in  the  terminal  elements. 
The  error  estimates  are  then  of  high  quality  in  the  entire  range  of  the 
accuracy . 

c)  Biquadratic  elements  lead  often  to  an  exponential  convergence  in 
the  entire  range  of  engineering  accuracy. 

d)  The  feedback  algorithm  should  be  based  on  a  properly  designed 
"short  passes"  which  anticipate  the  pattern  of  the  mesh  refinement. 
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